/**
   Copyright 2008 Kirill Ignatiev

    This file is part of kerr-wavefronts.

    kerr-wavefronts is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    kerr-wavefronts is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with kerr-wavefronts.  If not, see <http://www.gnu.org/licenses/>.
 */
#include <cstdio>
#include <cmath>
#include <vector>
#include <cassert>
#include <algorithm>
#include "solvegeods.h"
#include <gsl/gsl_odeiv.h>
using namespace std;



class geodesic {
public:
  double r;
  double phi;
  double t; // current time
  double data[3];

  gsl_odeiv_system sys;
  gsl_odeiv_control* ctl;
  gsl_odeiv_evolve* ev;
  gsl_odeiv_step* step;

  double start_x, start_y; // initial position

  geodesic(double r1, double phi1) : r(r1), phi(phi1), t(0) {
    sys.dimension = 3;
    ctl = gsl_odeiv_control_y_new(1e-6, 1e-6);
    ev = gsl_odeiv_evolve_alloc(3);
    step = gsl_odeiv_step_alloc(gsl_odeiv_step_rkf45, 3);
    start_x = r1*cos(phi1), start_y = r1*sin(phi1); // initial pos
  }

  virtual void evolve(double tmax) = 0;
  void evolve_ode(double tmax);
  virtual bool not_on_horizon() const = 0;
  double x() const { return r*cos(phi); }
  double y() const { return r*sin(phi); }
};

class schw_geodesic : public geodesic {
public:
  double l; // schw parameter

  schw_geodesic(double r0, double phi0, double rdot0, double phidot0);
  void evolve(double tmax);
  bool not_on_horizon() const { return r > 1+1e-3; }
};

class kerr_geodesic : public geodesic {
public:
  kerr_param param; // kerr parameter

  kerr_geodesic(const kerr_param& p, double r0, double phi0, double pr0);
  void evolve(double tmax);
  bool not_on_horizon() const { return r > param.rplus+1e-3; }
};



